
      module mo_schu

!!$      use ppgrid, only : pverp
!!$      use mo_grid, only : plevp
      use shr_kind_mod, only: r8 => shr_kind_r8

      implicit none

      integer, parameter :: ngast = 17
      integer, parameter :: tdim  = 501
      real(r8), parameter    :: t_del = 5./(tdim-1)
      real(r8), parameter    :: t_fac = (tdim-1)/5.

      integer :: ii, jj
      real(r8)    :: d_table(0:tdim,ngast-1)
      real(r8)    :: x_table(0:tdim,ngast-1)
      real(r8)    :: o2_table(tdim)
      real(r8), dimension(12,ngast-1) :: a, b

!-----------------------------------------------------------------------------
!  a(16,12)             coefficients for rj(m) (table 1 in kockarts 1994)
!  b(16,12)                              rj(o2)(table 2 in kockarts 1994)
!  rjm                  attenuation coefficients rj(m)
!  rjo2                 rj(o2)
!-----------------------------------------------------------------------------
      data ((a(jj,ii),jj=1,12),ii=1,ngast-1) / &
!a 57000-56500.5 cm-1
       1.13402e-01,1.00088e-20,3.48747e-01,2.76282e-20,3.47322e-01,1.01267e-19, &
       1.67351e-01,5.63588e-19,2.31433e-02,1.68267e-18,0.00000e+00,0.00000e+00, &
!a 56500-56000.5 cm-1
       2.55268e-03,1.64489e-21,1.85483e-01,2.03591e-21,2.60603e-01,4.62276e-21, &
       2.50337e-01,1.45106e-20,1.92340e-01,7.57381e-20,1.06363e-01,7.89634e-19, &
!a 56000-55500.5 cm-1
       4.21594e-03,8.46639e-22,8.91886e-02,1.12935e-21,2.21334e-01,1.67868e-21, &
       2.84446e-01,3.94782e-21,2.33442e-01,1.91554e-20,1.63433e-01,2.25346e-19, &
!a 55500-55000.5 cm-1
       3.93529e-03,6.79660e-22,4.46906e-02,9.00358e-22,1.33060e-01,1.55952e-21, &
       3.25506e-01,3.43763e-21,2.79405e-01,1.62086e-20,2.10316e-01,1.53883e-19, &
!a 55000-54500.5 cm-1
       2.60939e-03,2.33791e-22,2.08101e-02,3.21734e-22,1.67186e-01,5.77191e-22, &
       2.80694e-01,1.33362e-21,3.26867e-01,6.10533e-21,1.96539e-01,7.83142e-20, &
!a 54500-54000.5 cm-1
       9.33711e-03,1.32897e-22,3.63980e-02,1.78786e-22,1.46182e-01,3.38285e-22, &
       3.81762e-01,8.93773e-22,2.58549e-01,4.28115e-21,1.64773e-01,4.67537e-20, &
!a 54000-53500.5 cm-1
       9.51799e-03,1.00252e-22,3.26320e-02,1.33766e-22,1.45962e-01,2.64831e-22, &
       4.49823e-01,6.42879e-22,2.14207e-01,3.19594e-21,1.45616e-01,2.77182e-20, &
!a 53500-53000.5 cm-1
       7.87331e-03,3.38291e-23,6.91451e-02,4.77708e-23,1.29786e-01,8.30805e-23, &
       3.05103e-01,2.36167e-22,3.35007e-01,8.59109e-22,1.49766e-01,9.63516e-21, &
!a 53000-52500.5 cm-1
       6.92175e-02,1.56323e-23,1.44403e-01,3.03795e-23,2.94489e-01,1.13219e-22, &
       3.34773e-01,3.48121e-22,9.73632e-02,2.10693e-21,5.94308e-02,1.26195e-20, &
!a 52500-52000.5 cm-1
       1.47873e-01,8.62033e-24,3.15881e-01,3.51859e-23,4.08077e-01,1.90524e-22, &
       8.08029e-02,9.93062e-22,3.90399e-02,6.38738e-21,8.13330e-03,9.93644e-22, &
!a 52000-51500.5 cm-1
       1.50269e-01,1.02621e-23,2.39823e-01,3.48120e-23,3.56408e-01,1.69494e-22, &
       1.61277e-01,6.59294e-22,8.89713e-02,2.94571e-21,3.25063e-03,1.25548e-20, &
!a 51500-51000.5 cm-1
       2.55746e-01,8.49877e-24,2.94733e-01,2.06878e-23,2.86382e-01,9.30992e-23, &
       1.21011e-01,3.66239e-22,4.21105e-02,1.75700e-21,0.00000e+00,0.00000e+00, &
!a 51000-50500.5 cm-1
       5.40111e-01,7.36085e-24,2.93263e-01,2.46742e-23,1.63417e-01,1.37832e-22, &
       3.23781e-03,2.15052e-21,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
!a 50500-50000.5 cm-1
       8.18514e-01,7.17937e-24,1.82262e-01,4.17496e-23,0.00000e+00,0.00000e+00, &
       0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
!a 50000-49500.5 cm-1
       8.73680e-01,7.13444e-24,1.25583e-01,2.77819e-23,0.00000e+00,0.00000e+00, &
       0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
!a 49500-49000.5 cm-1
       3.32476e-04,7.00362e-24,9.89000e-01,6.99600e-24,0.00000e+00,0.00000e+00, &
       0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00 /


      data ((b(jj,ii),jj=1,12),ii=1,ngast-1) / &
!  57000-56500.5 cm-1
       1.07382e-21,9.95029e-21,7.19430e-21,2.48960e-20,2.53735e-20,7.54467e-20, &
       4.48987e-20,2.79981e-19,9.72535e-20,9.29745e-19,2.30892e-20,4.08009e-17, &
!  56500-56000.5 cm-1
       3.16903e-22,1.98251e-21,5.87326e-22,3.44057e-21,2.53094e-21,8.81484e-21, &
       8.82299e-21,4.17179e-20,2.64703e-20,2.43792e-19,8.73831e-20,1.46371e-18, &
!  56000-55500.5 cm-1
       1.64421e-23,9.26011e-22,2.73137e-22,1.33640e-21,9.79188e-22,2.99706e-21, &
       3.37768e-21,1.39438e-20,1.47898e-20,1.04322e-19,4.08014e-20,6.31023e-19, &
!  55500-55000.5 cm-1
       8.68729e-24,7.31056e-22,8.78313e-23,1.07173e-21,8.28170e-22,2.54986e-21, &
       2.57643e-21,9.42698e-21,9.92377e-21,5.21402e-20,3.34301e-20,2.91785e-19, &
!  55000-54500.5 cm-1
       1.20679e-24,2.44092e-22,2.64326e-23,4.03998e-22,2.53514e-22,8.53166e-22, &
       1.29834e-21,3.74482e-21,5.12103e-21,2.65798e-20,2.10948e-20,2.35315e-19, &
!  54500-54000.5 cm-1
       2.79656e-24,1.40820e-22,3.60824e-23,2.69510e-22,4.02850e-22,8.83735e-22, &
       1.77198e-21,6.60221e-21,9.60992e-21,8.13558e-20,4.95591e-21,1.22858e-17, &
!  54000-53500.5 cm-1
       2.36959e-24,1.07535e-22,2.83333e-23,2.16789e-22,3.35242e-22,6.42753e-22, &
       1.26395e-21,5.43183e-21,4.88083e-21,5.42670e-20,3.27481e-21,1.58264e-17, &
!  53500-53000.5 cm-1
       8.65018e-25,3.70310e-23,1.04351e-23,6.43574e-23,1.17431e-22,2.70904e-22, &
       4.88705e-22,1.65505e-21,2.19776e-21,2.71172e-20,2.65257e-21,2.13945e-17, &
!  53000-52500.5 cm-1
       9.63263e-25,1.54249e-23,4.78065e-24,2.97642e-23,6.40637e-23,1.46464e-22, &
       1.82634e-22,7.12786e-22,1.64805e-21,2.37376e-17,9.33059e-22,1.13741e-20, &
!  52500-52000.5 cm-1
       1.08414e-24,8.37560e-24,9.15550e-24,2.99295e-23,9.38405e-23,1.95845e-22, &
       2.84356e-22,3.39699e-21,1.94524e-22,2.72227e-19,1.18924e-21,3.20246e-17, &
!  52000-51500.5 cm-1
       1.52817e-24,1.01885e-23,1.22946e-23,4.16517e-23,9.01287e-23,2.34869e-22, &
       1.93510e-22,1.44956e-21,1.81051e-22,5.17773e-21,9.82059e-22,6.22768e-17, &
!  51500-51000.5 cm-1
       2.12813e-24,8.48035e-24,5.23338e-24,1.93052e-23,1.99464e-23,7.48997e-23, &
       4.96642e-22,6.15691e-17,4.47504e-23,2.76004e-22,8.26788e-23,1.65278e-21, &
!  51000-50500.5 cm-1
       3.81336e-24,7.32307e-24,5.60549e-24,2.04651e-23,3.36883e-22,6.15708e-17, &
       2.09877e-23,1.07474e-22,9.13562e-24,8.41252e-22,0.00000e+00,0.00000e+00, &
!  50500-50000.5 cm-1
       5.75373e-24,7.15986e-24,5.90031e-24,3.05375e-23,2.97196e-22,8.92000e-17, &
       8.55920e-24,1.66709e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
!  50000-49500.5 cm-1
       6.21281e-24,7.13108e-24,3.30780e-24,2.61196e-23,1.30783e-22,9.42550e-17, &
       2.69241e-24,1.46500e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, &
!  49500-49000.5 cm-1
       6.81118e-24,6.98767e-24,7.55667e-25,2.75124e-23,1.94044e-22,1.45019e-16, &
       1.92236e-24,3.73223e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00 /

      contains

      subroutine schu_inti
!-----------------------------------------------------------------------------
!	... initialize the tables
!-----------------------------------------------------------------------------

      implicit none

!-----------------------------------------------------------------------------
!	... local variables
!-----------------------------------------------------------------------------
      integer :: i, iw, k, j1, jp1, j
      real(r8)    :: col
      real(r8), dimension(6) :: a0, a1, b0, b1

      do iw = 1,ngast-1
	 x_table(0,iw) = sum( a(1:11:2,iw) )
	 d_table(0,iw) = sum( b(1:11:2,iw) )
         do k = 1,tdim
	    col = 22. + t_del*real(k-1)
	    o2_table(k) = col
	    col = 10.**col
	    a1(:) = a(2:12:2,iw)*col
	    b1(:) = b(2:12:2,iw)*col
	    where( a1(:) < 500. )
	       a0(:) = exp( -a1(:) )
	    elsewhere
	       a0(:) = 0.
	    endwhere
	    where( b1(:) < 500. )
	       b0(:) = exp( -b1(:) )
	    elsewhere
	       b0(:) = 0.
	    endwhere
	    x_table(k,iw) = dot_product( a(1:11:2,iw),a0(:) )
	    d_table(k,iw) = dot_product( b(1:11:2,iw),b0(:) )
         end do
      end do

      end subroutine schu_inti

      subroutine schu( o2col, secchi, dto2, xscho2 )
!-----------------------------------------------------------------------------
!   purpose:
!   calculate the equivalent absorption cross section of o2 in the sr bands.
!   the algorithm is based on:  g.kockarts, penetration of solar radiation
!   in the schumann-runge bands of molecular oxygen:  a robust approximation,
!   annales geophysicae, v12, n12, pp. 1207ff, dec 1994.  calculation is
!   done on the wavelength grid used by kockarts (1994).  final values do
!   include effects from the herzberg continuum.
!-----------------------------------------------------------------------------
!   parameters:
!   nz      - integer, number of specified altitude levels in the working (i)
!             grid
!   o2col   - real, slant overhead o2 column (molec/cc) at each specified (i)
!             altitude
!   dto2    - real, optical depth due to o2 absorption at each specified  (o)
!             vertical layer at each specified wavelength
!   xscho2  - real, molecular absorption cross section in sr bands at     (o)
!             each specified wavelength.  includes herzberg continuum
!-----------------------------------------------------------------------------

!!$      use mo_grid, only : plev, plevp
      use ppgrid, only : pver,pverp

      implicit none

!-----------------------------------------------------------------------------
!	... dummy arguments
!-----------------------------------------------------------------------------
      real(r8), intent(inout) :: dto2(pver,ngast-1)
      real(r8), intent(inout) :: xscho2(pverp,ngast-1)
      real(r8), intent(in)    :: o2col(pverp)
      real(r8), intent(in)    :: secchi(pverp)

!-----------------------------------------------------------------------------
!	... local variables
!-----------------------------------------------------------------------------
      integer :: i, iw, j, j1, jp1, k, ki, kp1
      integer :: index(pverp)
      integer :: minind(1), maxind(1)
      real(r8)    :: a0, a1, b0, b1
      real(r8), dimension(6)        :: ac, bc
      real(r8), dimension(pverp,6)  :: aa, bb
      real(r8), dimension(pverp)    :: rjm, rjo2
      real(r8), dimension(pverp)    :: rjmi, rjo2i, lo2col, dels


!-----------------------------------------------------------------------------
! 	... initialize r(m)
!-----------------------------------------------------------------------------
      rjm(1:pverp)  = 0._r8
      rjo2(1:pverp) = 0._r8

!-----------------------------------------------------------------------------
! 	... initialize the table interpolation
!-----------------------------------------------------------------------------
      where( o2col(:) /= 0 )
         lo2col(:) = log10( o2col(:) )
      endwhere
      do ki = 1,pverp
	 if( o2col(ki) /= 0. ) then
	    if( lo2col(ki) <= o2_table(1) ) then
	       dels(ki)  = 0._r8
	       index(ki) = 1
	    else if( lo2col(ki) >= o2_table(tdim) ) then
	       dels(ki)  = 1._r8
	       index(ki) = tdim-1
	    else
               do k = 2,tdim
	          if( lo2col(ki) <= o2_table(k) ) then
	             dels(ki)  = t_fac*(lo2col(ki) - o2_table(k-1))
	             index(ki) = k-1
		     exit
		  end if
               end do
	    end if
	 else
	    index(ki) = 0
	    dels(ki)  = 0._r8
	 end if
      end do

!-----------------------------------------------------------------------------
! 	... calculate sum of exponentials (eqs 7 and 8 of kockarts 1994)
!-----------------------------------------------------------------------------
      do iw = 1,ngast-1
         do k = 1,pverp
	    ki = index(k)
	    rjm(k)  = x_table(ki,iw) + dels(k)*(x_table(ki+1,iw) - x_table(ki,iw))
	    rjo2(k) = d_table(ki,iw) + dels(k)*(d_table(ki+1,iw) - d_table(ki,iw))
	 end do

         do k = 1,pver
            if( rjm(k) > 1.e-100_r8 ) then
               kp1 = k + 1
               if( rjm(kp1) > 0. ) then
                  dto2(k,iw) = log( rjm(kp1) ) / secchi(kp1) - log( rjm(k) ) * secchi(k)
               else
                  dto2(k,iw) = 1000._r8
               end if
            else
               dto2(k,iw)   = 1000._r8
            end if
         end do
         do k = 1,pverp
            if( rjm(k) > 1.e-100_r8 ) then
               if( rjo2(k) > 1.e-100_r8 ) then
                  xscho2(k,iw) = rjo2(k)/rjm(k)
               else
                  xscho2(k,iw) = 0._r8
               end if
            else
               xscho2(k,iw) = 0._r8
            end if
         end do
      end do

      end subroutine schu

      end module mo_schu
